0:50 Inverse Problems in Science and Engineering DT'SPIE'2012'8'20 



Inverse Problems in Science and Engineering 
Vol. 00, No. 00, August 2012, 1-23 



RESEARCH ARTICLE 

Optical imaging of phantoms from real data by an approximately 
globally convergent inverse algorithm 

Jianzhong Su^, Michael V. Klibanov'', Yueming Liu^, Zhijin Lin'^, Natee Pantong'^ and 

Hanli Liu'^ 

'^Department of Mathematics, University of Texas at Arlington, Arlington, TX 76019, 
USA; ^Department of Mathematics and Statistics, University of North Carolina at 
Charlotte, Charlotte, NC 28223, USA'^ Department of Bioengineering, University of Texas 
at Arlington, Arlington, TX 76019, US A'^ Department of Mathematics and Computer 
Science, Royal Thai Air Force Academy, Bangkok, Thailand 
(Received 00 Month 200x; in final form 00 Month 200x) 

A numerical method for an inverse problem for an elliptic equation with the running source 
at multiple positions is presented. This algorithm does not rely on a good first guess for the 
solution. The so-called "approximate global convergence" property of this method is shown 
here. The performance of the algorithm is verified on real data for Diffusion Optical Tomog- 
raphy. Direct applications are in near-infrared laser imaging technology for stroke detection 
in brains of small animals. 



Keywords: Approximate Global Convergence Property; Inverse Problem; Diffusion Optical 
Tomography; Real Data 

AMS Subject Classification: 65B21; 65D10; 65F10 



1. Introduction 

We consider a Coefficient Inverse Problem (CIP) for a partial differential equation 
(PDE) - the diffusion model with the unknown potential. The boundary data for 
this CIP, which model measurements, are originated by a point source running 
along a part of a straight line. This PDE governs light propagation in a diffusive 
medium, such as, e.g. biological tissue, smog, etc.. Thus, our CIP is one of problems 
of Diffusion Optical Tomography (DOT). We are interested in applications of DOT 
to the detection of stroke in small animals using measurements of near infrared light 
originated by lasers. Hence, the above point source is the light source in our case. 
The motivation of imaging of small animals comes from the idea that it might be a 
model case for the future stroke detection in humans, at least in brains of neonatals, 
via DOT. We apply our numerical method to a set of real data for a phantom 
medium modeling the mouse's brain. Although this algorithm was developed in 



earlier publications [17|, I24l426l| of this group, its experimental verification is new 
and is the main contribution of this paper. 

As to our numerical method, we introduce a new concept of the "approximate 
global convergence" property. In the previous publication [17| of this group about 
this method the approximate global convergence property was established in the 



continues case. Compared with 17|, the main new analytical result here is that 



we establish this property for the more realistic discrete case. In the convergence 
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analysis of [17] the Schauder theorem [14.] was apphed for C^"*""— solutions of cer- 
tain elliptic equations arising in our method. Now, however, since we consider the 
discrete case, we use the Lax-Milgram theorem. Here and below 0^+°^ are Holder 
spaces fl2], where /c > is an integer and a G (0, 1). 

CIPs are both nonlinear and ill-posed. These two factors cause very serious chal- 
lenges in their numerical treatments. Indeed, corresponding least squares Tikhonov 
regularization functionals usually suffer from multiple local minima and ravines. 
As a result, conventional numerical methods for CIPs are locally convergent ones, 
see, e.g. j3] and references cited there. To have a guaranteed convergence to a true 
solution, a locally convergent algorithm should start from a point which is located 
in a small neighborhood of this solution. However, it is a rare case in applications 
when such a point is known. The main reason why our method avoids local minima 
is that it uses the structure of the underlying PDE operator and does not use a 
least squares functional. 

Therefore, it is important to develop such numerical methods for CIPs, which 
would have a rigorous guarantee of providing a point in a small neighborhood of 
that exact solution without any a priori knowledge of this neighborhood. It is 
well known, however, that the goal of the development of such algorithms is a 
substantially challenging one. Hence, it is unlikely that such numerical methods 
would not rely on some approximations. This is the reason why the notion of the 
approximate global convergence property was introduced in the recent work [l^ : 
also see section 1.1.2 of the book [9\ and subsection 4.1 below. The verification 
of our approximately globally convergent numerical method on computationally 



simulated data was done in [17|, I24l426f] . In this paper we make the next step: 



we verify this method on real data. Regardless on some approximations we have 
here, the main point is that our numerical method does not rely on any a priori 
knowledge of a small neighborhood of the exact solution. 
In 0-0; 3" 13] a similar numerical method for CIPs for a hyperbolic PDE was 



developed. These results were summarized in the boo k M . In particular, it was 
demonstrated numerically in Test 5 of IQ], section 8.4 of [ig ]. on page 185 of [1] and 
in section 5.8.4 of [§] that the approximately globally convergent method of 0- 
[s, 18-20] outperforms such locally convergent algorithms as quasi-Newton method 
and gradient method. The main difference of the technique of these publications 
with the one of the current paper is in the truncation of a certain integral. In 



1814201] that integral was truncated at a high value s > of the parameter 
s > of the Laplace transform of the solution of the underlying hyperbolic PDE. 
Indeed, in the hyperbolic case the truncated residual of that integral, which we call 
the "tail function", is 0(l/s), as s — )• oo, i.e. this residual is automatically small 
in this case. Being different from the hyperbolic case, in the elliptic PDE with 
the source running along a straight line, s > represents the distance between 
the source and the origin. Because of this, the tail function is not automatically 
small in our case. Thus, a special effort to ensure this smallness was undertaken in 
17], and is repeated in the current publication. It was shown numerically in [l3] 



that this special effort indeed improves the quality of the reconstruction, compare 



Figures 2b and 2c in |l7l ]. 



Our real data were collected at the boundary of a 2-D cross-section of the 3-D 
domain of interest, see section 6. Hence, we have imaged this cross-section only 
and have ignored the dependence on the third variable. In addition, our theory 
requires the use of many sources placed along a straight line. However, limitations 
of our experimental device allow us to use only three sources on this line. We 
believe that the accuracy of our results (Table 8.2) confirms the validity of both 
these assumptions as well as manifests a well known fact that there are always 
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some discrepancies between the analysis and computational studies of numerical 
methods. 

Because of the above mentioned substantial challenges, the topic of the devel- 
opment of non-locally convergent numerical methods for CIPs is currently in its 
infancy. As to such methods for CIPs for elliptic PDEs, we refer to, e.g. publica- 
tions 0,[il|,[il,[2i|-[2i and references cited there. These publications are concerned 



with the Dirichlet-to-Neumann map (DN). We are not using the DN here. 

In a general setting, an ill-posed problem is a problem of solving an equation 
F (x) = y with a compact operator F. Here x G Bi, y £ B2, where Bi and B2 
are certain Banach spaces. Naturally, the operator F varies from one problem to 
another one. It is well known that the existence theorem for such an equation is very 
tough to prove, since the range of any compact operator is 'too narrow' in certain 
sense, see, e.g. Theorem 1.2 in 0]. Therefore, by one of fundamental concepts of the 
theory of Ill-Posed Problems, one should assume the existence of an exact solution 
of such a problem for the case of an "ideal" noiseless data @, 0, ll^j- Although this 
solution is never known in practice and is never achieved in practice (because of the 
noise in the real data), the regularization theory says that one needs to construct 
a good approximation for it. We assume that our CIP has an exact solution for 
noiseless data and also assume that this solution is unique. 

The rest of this paper is arranged as follows. In section 2 we pose both forward 
and inverse problems and study some properties of the solution of the forward 
problem. In section 3 we present our numerical method. In section 4 we conduct 
the convergence analysis. In section 5 we discuss the numerical implementation of 
our method. In section 6 we describe the experiment. In section 7 we outline our 
procedure of processing of real data. In section 8 we present reconstruction results. 
We briefly summarize results in section 9. 



2. Statement of the Problem 
2.1. The Inverse Problem 

Below X = (x, z) S and C is a convex bounded domain. The boundary of 
this domain dQ E in our analysis. In numerical studies dO. is piecewise smooth. 
Consider the following elliptic equation in with the solution vanishing at infinity, 

Au — a (x) u = —6 (x — xq) , x, xq G R^, (1) 
lim n(x,xo) = 0. (2) 

|x|— >CX) 

Inverse Problem. Let k = const. > 0. Suppose that in (Cp the coefficient a (x) 
satisfies the following conditions 

aeC^ (R2) , a (x) > and a (x) = k^ for x G R^\n. (3) 

Let L C (R^\r2) be a straight line and T C L be an unbounded and connected 
subset of L. Determine the function a (x) inside of the domain assuming that 
the constant k is given and also that the following function (/?(x, xq) is given 

U (X, Xq) = if (X, Xq) , V (x, Xq) G 5il X P. (4) 

We are unaware about a uniqueness theorem for this inverse problem. Neverthe- 
less, because of the above application, it is reasonable to we develop a numerical 
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method. Thus, we assume that u niq ueness for this problem holds. Our numerical 



studies of both the past jlTl , l24l426l| and the current publication indicate that a 
certain uniqueness theorem might be established. 

We assume that sources {xq} are located outside of the domain of interest Q. 
because this is the case of our measurements and because we do not want to work 
with singularities in our numerical method. The CIP ([l])-([4]) has an application 
in imaging using light propagation in a diffuse medium, such as biological tissues. 
Since the modulated frequency equals zero in our case, then this is the so-called 
continuous-wave (CW) light. The coefficient a (x) = 3 {fJ-'sfJ-a) (x) , where fi'^ (x) is 
the reduced scattering coefficient and fia (x) is the absorption coefficient of the 
medium [1, 0] • In the case of our particular interest in stroke detections in brains 
of small animals, the area of an early stroke can be modeled as a small sharp 
inclusion in an otherwise slowly fluctuating background. Usually the inclusion/ 
background contrast aind/cLb > 2. Therefore our focus is on the reconstruction of 
both locations of sharp small inclusions and the values of the coefficient a (x) inside 
of them, rather than on the reconstruction of slow changing background functions. 



2.2. Some Properties of the Solution of the Forward Problem ^\), ([![) 

2.2.1. Existence and uniqueness 

First, we state the existence and uniqueness of the solution of the forward prob- 
lem ([1]), ([2]). For brevity we consider only the case xq ^ 17, since this is the case 
of our Inverse Problem. Let Kp (z) ,z G [R,p > be the Macdonald function. It is 
well known [1,] that for y G K 



Kp{y) = ^e''y[l + 0[-] ),y^oo. (5) 



2y 

Theorem 2.1 : Let $7 C be the above bounded domain. Assume that the 
coefficient a (x) satisfies conditions ([3]) . Then for each source position xq G R^\r2 
there exists unique solution u(x, Xq) of the problem ([T|), ([2|) such that 

^(x,Xo) = 7^-fCo(fc|x-xo|)-Fu(x,Xo) :=uo(x-xo)-FS(x,Xo), (6) 

ZTT 

where the function uq is the fundamental solution of equation ([1]) with a (x) = 
fc^, the function u satisfies Q, u G i/^ (R^) and u G 0^+" (R^) . In addition, 
u (x, Xq) > 0, Vx G Jl. 



A similar result for the 3-D case was proven in [9|, see Theorem 2.7.2 in this 
reference. Hence, we leave out the proof of Theorem 2.1 for brevity. 

2.2.2. The asymptotic behavior at |xo| — )• oo 

It follows from ([5]) and ^ that the asymptotic behavior of the function 
Uo (x - Xo) is 

■Uo (x - Xo) = wq (|xo|) (^ + (^1^^^ ) |xo| OO, 

g-fc|xo| 

^o(|xo|) = pT^r 
2y^27r|xo| 
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Denote 6 (x) = a (x) - P. Then by ([3]) 6 (x) = for x G R'^\Q. Let Mi be a 
positive constant. Denote 

B{Mi) = |6 G (R^) : < Mi,6(x) > 0,6(x) = for x G R^^^J)} . 

Also, let the function p^o (x) satisfies conditions 

Poo (x) G C72+° (|x| <R),yR> 0,poo G (r2) . (7) 

and be the solution of the following problem 

Apoo-A;^Poo-&(x)poo = &(x),xG R^ (8) 
hm Poo (x) = 0. (9) 

\x\^oo 

The uniqueness and existence of the solution of the problem ©-(IH]) are similar to 
these of Theorem 12. 1[ 

Lemma 2.2: 1 +Poo > 0. 

Proof: Let p = 1 + Poo • Then 

Ap — k'^p — b{x)p = —k'^, lim p(x) = l. (10) 

|x|— >oo 

Consider a sufficiently large number R > such that p (x) > 1/2 for x G {|x| > R} . 
Then the maximum principle [14] applied to equation (jlOp for x G {|x| < R} shows 
that p(x) > in {|x| < i?}. □ 



Lemma 2.3: JX3/ Let the function b G B(Mi). Then there exists a constant 
M2 (Ml, n) > such that 

||lnu(x,Xo) - Inwo (|xo|) - ln(l +poo {^))\\c^m) < \ — \ , 

^ ' Xq 



Xo G {|xo| > 1} n {R^\n) ,x G 17. 



3. The Numerical Method 

Since we can put the origin on the straight line L, if necessary, then without any loss 
of the generality we can set s := |xo|, assuming that only the parameter s changes 
when the source xq runs along T C L. Denote u (x, s) := u (x, xq) , x G il, xq G L. 
Since L n O = and the point source xq G T, then xq ^ Since by Theorem 
I2.1l u (x, s) > 0, Vx G ri, then we can consider the function w (x, s) =lnu (x, s) for 
X G O. We obtain from ^ and (g]) 

Aw + \Vw\'^ = a{x) inn, (11) 
w (x, s) = ifi (x, s) , V (x, s) G 90 X [s, s] , (12) 
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where 991 = In and s, s are two positive numbers, which should be chosen in 
numerical experiments. 



3.1. The integral differential equation 

We now "ehminate" the coefficient a (x) from equation (jlip via the differentiation 
with respect to s. However, to make sure that the resulting the so-called "tail 
function" is small, we use the above mentioned (Introduction) special effort of the 
paper [13]. Namely, we divide pT]) by s^^p > 0. In principle, any number p > can 
be used. But since in computations we took p = 2, both here and in [13], then we 
use below only p = 2, for the sake of definiteness. Denote w (x, s) = w (x, s) / s^. 
By pT]) equation for the function w (x, s) is 



Aw + s''\Vw\= (13) 
Next, let g (x, s) := S^u) (x, s) = ds Inn (x, s)) , for s G [s, s]. Then 

s 

■w {x, s) = — J g (x, r) dr + T (x) , X G il, s G [s, s] , (14) 

s 

where T (x) is the so-called "tail function" . The exact expression for this function is 
of course T (x) = w (x, s). However, since the function w (x, s) is unknown, we will 



use an approximation for the tail function, see subsection 5.2 as well as [17l]. By the 
Tikhonov concept for ill-posed problems , one should have some a priori 

information about the solution of an ill-posed problem. Thus, we can assume the 
knowledge of a constant Mi > such that the function a (x) — /c^ = b (x) G B (Mi). 
Hence, it follows from Lemma 12.31 that 

T(x,^) = i^^ + i^^ii±|^ + ^, xGf^, V^>1, (15) 
s s s 

\\g {x,s)\y(Ji) <M2(Mi,0), Vs> 1, V6gS(Mi), 

where the number M2 (Mi, fi) is independent on s. Differentiating ()13p with respect 
to s, we obtain the following nonlinear integral differential equation for the function 



q m 



Ag - ^ y Aq (x, r) dr - 2s^Vq j Vq (x, r) dr (16) 



+4sl-j Vq (x, r) + VT j + 2s^yTVq = -^AT. 



In addition, (fT^ implies that 



q (x, s) = 'il) (x, s) , V (x, s) edVtx [s, s] , (17) 
V'(x,s) =5, (s-2lnv9(x,s)) . (18) 
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If we approximate well both functions q and T together with their derivatives up 
to the second order, then we can also approximate well the target coefficient a (x) 
via backwards calculations, see (I4ip . Therefore, the main questions now is: How to 
approximate well both functions q and T using Ij^l^^-^T^)? 



3.2. Layer stripping with respect to the source position 

In this subsection we present a layer stripping procedure with respect to s for 
approximating the function q, assuming that the function T is known. Usually the 
layer stripping procedure is applied with respect to a spatial variable. However, 
sometimes the presence of a differential operator with respect to this variable in 
the underlying PDE results in the instability, since computing derivatives is an 
unstable procedure. The reason why our layer stripping procedure is stable is that 
we do not have a differential operator with respect to s in our PDE. 

We approximate the function q (x, s) as a piecewise constant function with re- 
spect to the source position s. We assume that there exists a partition s = sn< 
sm-i < . . . < Si < sq = s, Sj-i — Si = h oi the interval [s,s] with a sufficiently 
small step size h such that 

q (x, s) = qn (x) for s £ [s„, s„_i) , n > 1; qo := 0. (19) 

Hence, 

^ n-1 

/ q (x, s) ds = {sn-i - s) qn (x) + /i ^ qj (x) . 

Let ipn (x) be the average of the function ip (x) over the interval {sn, Sn-i)- Then 
we approximate the boundary condition (I17p as a piecewise constant function with 
respect to s, 

g„(x) = ?/;„(x),xe517. (20) 

Using integrate equation (fTUj) with respect to s G [sn,s„_i). We obtain for 
n > 1 

Aqn + A2,n ( h J2 ^Qj - VT ) Vqn - Ai^n {yqnf = 

n— 1 / n— 1 \ 

A^^nh E + A4.„ U E - VT - ^3,nAr, 
i=i V i=o / 



where A^^n, k = 1, 4 are certain coefficients, see [l7l] for exact formulas for them. 
Let 

s > 2,/i e (0, 1). (22) 
Then one can prove that |^j,n| < 8s^,i = 2,3,4, 

l^i.nl < 2s^/i. (23) 
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Assuming that h is sufficiently small and using ()23p . we assume below that 

^i,n(V<7n)' :=0 in ([21]). (24) 

It is convenient for our convergence analysis to formulate the Dirichlet boundary 
value problem (j20p . (j2ip in the weak form. Denote HQ{yt) = {li S H^{yt) : u \dn= 
0}. Assume that there exists functions such that 

G i^) and ^nlgn = Vn (x) , n G [1, N] . (25) 

Consider the function Pn = Qn — G Hq (0) . Then we obtain an obvious analog 
of equation ([2T]) for the function p„. Multiply both sides of the latter equation by 
an arbitrary function 7] G Hq ($7) and integrate over Q using integration by parts 
and (El). We obtain 



n-l 



J VpnVrjdx - A2^n J I ^ ^ Vgj - VT J Vp„ • r]dx 

ri}dx (26) 




i=o 

n-l 



j=0 



/n—1 
h J2 VgjV?7dx, Pn G (0) , Vry G //q^ (0) 

,o 3=0 



/„ (x) = ^4,n I /i X] "79^ - VT I - ^3,nAT. 

3=0 

Hence, the function g„ G (il) is a weak solution of the problem (pOl) . (i2T]l 
if and only if the function pn G Hq (Q.) satisfies the integral identity (f26l) . The 
question about existence and uniqueness of the weak solution of the problem (f26]l 
is addressed in Theorem 14.21 

We now describe our algorithm of sequential solutions of boundary value prob- 
lems (I20p . (I2ip for n = 1, . . . , N , assuming that an approximation T (x) for the tail 
function is found (see subsection 5.2 for the latter). Recall that qo = 0. Hence, we 
have: 

Step n G [1,-/V]- Suppose that functions qi, ■ ■ ■ ,qn-i are computed. On this step 
we find the weak solution of the Dirichlet boundary problem (I20p . (I2ip for the 
function g„ via the FEM with triangular finite elements. 

Step + 1. After functions qi,. . . ,qN are computed, the function a (x) is re- 
constructed using backwards calculations at n := N, i.e., for the lowest value of 
s := S]y = s as described in subsection 5.1. 
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4. Convergence Analysis 

4.1. Approximate global convergence 



The central question which was addressed in above cited pubUcations [17|, l24l426l| 
was of constructing such a numerical method for the CIP ((l])-([4]), which would 
simultaneously satisfy the following three conditions: 

1. This method should deliver a good approximation for the exact solution of 
this CIP without an a priori knowledge of a small neighborhood of this solution. 

2. A theorem should be proven which would guarantee of obtaining such an 
approximation. 

3. A good numerical performance of this technique should be demonstrated on 
computationally simulated data and, optionally, on real data. 

We use the word "optionally" , because it is usually hard and expensive to actu- 
ally collect real data. It is enormously challenging to construct such a numerical 
method. Therefore, as it is often done in mathematical modeling, conceptually, our 
approach consists of the following six steps: 

Step 1. A reasonable approximate mathematical model is proposed. The accu- 
racy of this model cannot be rigorously estimated. 

Step 2. A numerical method is developed, which works within the framework 
of this model. 

Step 3. A theorem is proven, which guarantees that, within the framework of this 
model, the numerical method of Step 2 indeed reaches a sufficiently small neighbor- 
hood of the exact solution. Naturally, the smallness of this neighborhood depends 
on the level of the error, both in the data and in some additional approximations 
is sufficiently small. 

Step 4. The numerical method of Step 2 is tested on computationally simulated 
data. 

Step 5. The numerical method of Step 2 is tested on real data. 
Step 6. Finally, if results of Steps 4 and 5 are good ones, then that approximate 
mathematical model is proclaimed as a valid one. 
Thus, results of the current paper (Step 5) in combination with the previous ones 



of [17|, I24l426f| lead to the positive conclusion of Step 6. Step 6 is logical, because 
its condition is that the resulting numerical method is proved to be effective. It 
is sufficient to achieve that small neighborhood of the exact solution after a fi- 
nite (rather than infinite) number of iterations. We refer to page 157 of the book 
(l3l | where it is stated that the number of iterations can be regarded as a regular- 
ization parameter sometimes for an ill-posed problem. Still, in our computations 
both for simulated and real data, classical convergence in the Cauchy sense was 
achieved. These consideration lead to the following definition of the approximate 
global convergence property. 

Definition 4.1: (approximate global convergence) [^,[3]. Consider a nonlinear 
ill-posed problem P. Suppose that this problem has a unique solution x* £ B for the 
noiseless data y* , where S is a Banach space with the norm ||-||^ . We call x* "exact 
solution" or "correct solution" . Suppose that a certain approximate mathematical 
model M is proposed to solve the problem P numerically. Assume that, within 
the framework of the model M, this problem has unique exact solution x^. Also, 

let one of assumptions of the model M be that x\.j^ = x*. Consider an iterative 
numerical method for solving the problem P. Suppose that this method produces 
a sequence of points {xn}n=i C B, where N G [l,oo) . Let the number 9 G (0, 1) . 
We call this numerical method approximately globally convergent of the level 9, or 
shortly globally convergent, if, within the framework of the approximate model M, 
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a theorem is proven, which guarantees that, without any a priori knowledge of a 
sufficiently small neighborhood of x*, there exists a number N G [1, N) such that 

\\xn-x*\\j^<e,'ine\N,N]. (27) 

Suppose that iterations are stopped at a certain number k G [A^, A^] . Then the 
point Xfc is denoted as := Xgiob and is called "the approximate solution resulting 
from this method" . 

With reference to the notion of the approximate mathematical model M, as well 
as to the above Steps 1-6, it is worthy to mention here that one of the keys to the 
successful numerical implementation [2] of the non-local reconstruction algorithm 
of [12] was the use of a certain approximate mathematical model. The same is 
true for the two dimensional analog of the Gel'fand-Levitan-Krein equation being 
applied in [l^ to an inverse hyperbolic problem. Thus, it seems that whenever 
one is trying to construct an efficient non-locally convergent numerical method 
for a truly complicated nonlinear ill-posed problem, it is close to the necessity to 
introduce some approximations which cannot be rigorously justified and then work 
within the resulting approximate mathematical model then. Analogously, although 
the Huygens-Fresnel theory of optics is not rigorously supported by the Maxwell's 
system (see section 8.1 of the classical textbook [iq|), the "diffraction part" of the 
entire modern optical industry is based on the Huygens-Fresnel optics. 



4.2. Exact solution 

By one of statements of Introduction as well as by Definition 4.1, we should assume 
the existence and uniqueness of an "ideal" exact solution a* (x) of our inverse 
problem for an "ideal" noiseless exact data (p* (x, Xg) in Next, in accordance 
with the regularization theory, one should assume the presence of an error in the 
data of a small level 7 > and construct an approximate solution for this case. 



Since the exact solution was defined in [17| , we outline it only briefly here for the 



convenience of the reader. Let the function a* (x) satisfying conditions ([3]) be the 
exact solution of our inverse problem for the noiseless data ^p* (x, xq) in We 
assume that a* (x) is unique. Let the function u* (x, s) be the same as the function 
u (x,xo) in Theorem 12. H but for the case a (x) := a* (x). Denote 

w* (x, s) = s^^ \nu* (x, s) , q* (x, s) = dgW* (x, s) , 
T* (x) =w* (x,s). 



The function q* satisfies the analogue of equation (jl6p with the boundary condition 

m 

q* (x, s) = ij* (x, s) , V (x, s) ednx [s, s] , (28) 

where by ([H]) V* (x, s) = ds (s"^ Inu* (x, s)) for (x, s) G 9Qx [s, s]. We call q* (x, s) 
the exact solution. It should be noted that our real data f (x,s) are naturally given 
with a random noise. It is well known that the differentiation of the noisy data is 
an ill-posed problem [^,[23]. However, because of some limitations of our device, we 
work only with three values of the parameter s in our real data (subsection 5.2). 
Hence, we simply calculate the s— derivative via the finite difference. Our numerical 
experience shows that this does not lead to a degradation of our reconstruction 
results. 
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4.3. Our approximate mathematical model 

In this model we basically assume that the upper bound for the tail function can 
become sufficiently small independently on s. In addition, since we calculate the 
tail function T (x, s) separately from the function q, we assume that the function 
T (x, s) is given (but not the exact tail function T* (x, s)). Although the smallness of 
T is supported by ([TS]) . the independence of that upper bound of this function from 
s does not follow from (|15p . Still this is one of two assumptions of our approximate 
mathematical model. Following the concept of Steps 1-6 of subsection 4.1, we now 
verify this model on real data. 

More precisely, our approximate mathematical model M consists of the following 
two assumptions 

Assumptions: 

1. We assume that the number s > 1 is sufficiently large and fixed. Also, the tail 
function T (x, s) is given, and tail functions T* (x, s) , T (x, s) have the forms 

T* (x, s) = h^^^ + r* (x, s) , X G f), r* e (O) , (29) 

r(x,s) = i^;^^^^+r(x,s), ^en, reC^n), (30) 

\\r*\\c2(n) M\c-{n) < (31) 

where ^ G (0, 1) is a small parameter. Furthermore, the parameter ^ is independent 
on s. 

2. As the unknown vector x & B in Definition 4.1, we choose the vector function 
g(x) = {qi, ...,qN) (x) in (fT9]) . (f2T]) . In this case we will have only one iteration in 
Definition 4.1, i.e. iV = iV = 1 in ([27]) . 

It is for the sake of our convergence analysis that we choose here the vector func- 
tion q (x) rather than the coefficient a (x) as our unknown function. Indeed, finding 
a (x) in the weak form (j4ip . (j42p would require some interpolation estimates of the 
differences between functions and their interpolations via the finite dimensional 
subspace Gm (see subsection 4.4 for p* , Gm)- The latter has an "underwater rock" 
in terms of the convergence analysis. 

By Definition 4.1 we need to prove uniqueness of the function q* (x, s) for x € 
Q,s £ [s, s] . Assuming that the exact tail function T* (x,s) satisfying (p9]) . ([3T]) is 
given and the interval [s, s] is sufficiently small, this can be done easily, using the 
fact that q* (x, s) is the solution of an obvious analog of the problem (fT6]l . (fTTl) . 
This is because Volterra-like integrals are involved in ()16p . 

4.4. Approximate global convergence theorem 

Let 7 be the sum of the following errors: the level of the error in the function 
s) (by (jlSp this function is generated by the measurement data), errors in our 
finite element approximations of qn, the magnitude of the C^(r2) norm of the tail 
function as well as the length of the interval (s, s). 

Consider a finite dimensional subspace Gm C Hq (Q,) with dim Gm = m. A 
realistic example of Gm is the subspace, generated by triangular finite elements, 
which are piecewise linear functions. We assume that fx,fz £ L^o (f^) ,V/ G Gm- 
Since all norms in the subspace Gm are equivalent, then there exists a constant 
Gm = Gm (Gm) such that 

||V/||^^(^) < Gm ||V/||^^(^) , V/ G Gm, Gm > 2. (32) 
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1/2 



/ 2 2 \ ' 

Here and below ||V/||^^(j^) := [\\fx\\L^(n) + WfyWh^io.)) ' ^^^^ s^™^ ^^^^ 
L2 {Q) norm. Theorem 14.21 works with the weak formulation (j26p . In the course of 
the proof of this theorem we need to estimate products VgjVp„. While this was done 
for C^"*"" {yt) solutions in 13| using Schauder theorem, in the weak formulation we 



need to assume that functions 

qn - ^'n := Pn G Gm- (33) 

We introduce functions g*,^'*, which are analogs of functions qn,^n for the case 
a := a*. Let p* := - . We assume that 

q:,K e m , l|Vg:ile(TT) < C*,n G [1,N] , (34) 



max 

sS[s„,s„. 



(x,s) - g; (x)||^,(^)) < C*h,nG [1,7V] , (35) 



Ai,n {^qnf := 0, (36) 



p*^eGm,\\Vp*Jc(n)<'^C*,ne[l,N], (37) 



II VM^n (x) - (x)ll^^(f,) <C*{a + h),ne [1, iV] , (38) 

where the number C* > 1 is given and o" > is a small parameter characterizing 
the level of the error in the data ip (x, s). The connection between a and the error 
in ip{x,s) is clear from comparison of (pO|) . (p5|) with p8]) . Condition (j36|) is a 
direct analog of (|2l|) . Each function is the weak solution of an analog of the 
problem (p6|) . 



Theorem 4.2 : Let C be a convex bounded domain with the boundary dQ G 
C^. Assume that Assumptions 1,2 of subsection 4.3 hold. In addition, assume that 
conditions (f22]l . (fM|l . (p3]l - ([38]l hold. Let f3 = s — s. Consider the error parameter 
7 = /3 + /i + (T + ^, where is the small parameter defined in l^29\) -[31\ ) . which is 
independent on s (Assumption 1). Let Cm > 1 be the constant defined in ()32p . 
Then there exists a constant B = B (J7) > such that if 

then for each n G [1, A^] there exists unique solution p„ G Gm of the problem ((26]) 
and for functions qn= Pn + the following estimate holds 

l|V(Z„- Vg:||^^(f^) <0, (40) 

ly/iere 6* = {BC*s'^) 7 G (0, 1) . T/jms, (gU]) imphes that our numerical method has 
the approximate global convergence property of the level 9 (Definition 4-1)- 



Remarks 4.1. 
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1. We omit the proof of this theorem, because it is similar with the proof of a 
similar theorem in for the C" space. The only technical challenge in theorem 
14.21 is to obtain analogous arguments in the finite dimensional space Gm- 

2. The smallness of the parameter /3 in (j39p is a natural requirement since the 
original equation (I16p contains Volterra-like integrals in nonlinear terms. It is well 
known from the standard ODE course that the existence of a Volterra-like nonlinear 
integral equation of the second kind can be proven only on a small interval. 



5. Numerical Implementation 



Since this implementation was described in detail in [17[, we outline it only briefly 
here for the convenience of the reader. As to the functions qn, we have sequentially 
calculated them via the FEM solving Dirichlet boundary value problems (j20p . (j2ip . 
As it was mentioned in the end of subsection 3.2, we have used standard triangular 
finite elements. Two important questions which are discussed in this section are 
about approximating the function a (x) and the tail function T (x) . 



5.1. Approximation of the function a (x) 

We reconstruct the target coefficient a (x) via backwards calculations as follows. 
First, we reconstruct the function w {x, sn) from qn,n = 1,2, N and tail T. We 
have tt(x. Sat) = exp [s'^w {x, sn)]- Next, since in ([1]) the source xq ^ we use 
equation ([1]) in the weak form as 

— J Vu'Vr]pdx = J auijpdx, (41) 



where the test function rjp (x) ,p £ [1, P] is a quadratic finite element of a compu- 
tational mesh with the boundary condition r]p (x) =0. The number P is finite 
and depends on the mesh we choose. Equalities (I4ip lead to a linear algebraic sys- 
tem which we solve. Since this formulation is complex but standard, it is omitted 



here. Interested readers can see |l2l ]. Finally, we let 



a (x) = max (a (x) , /c^) . (42) 



5.2. Construction of the tail function 

The above construction of functions qn depends on the tail function T. In this 
subsection we state briefly our procedure of approximating the tail function, see 



171 . |24| . 1251 ] for details. This procedure consists of two stages. First, we find a first 
guess for the tail using the asymptotic behavior of the solution of the problem ([T]), 
([5]) as |xo| — )• oo (Lemma 12. 3p . as well as boundary measurements. On the second 
stage we refine the tail. 
We display in Fig. [1] the boundary data collection scheme in our experiment. We 



have used six locations of the light sources, see Fig. 1(a) Sources number 1,2 and 
3 are the ones which model the source xq running along the straight line L, see (j3|). 
The distance between these sources was six (6) millimeters. Each light source means 
drilling a small hole in the phantom to fix the source position. It was impossible to 
place more sources in a phantom of this size that mimics actual mouse head. So, 
since the data for the functions qn are obtained via the differentiation with respect 
to the source position, we have used only two functions qi,q2- Fortunately, the 
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(a) (b) 

Figure 1. The schematic diagrams of inverse problem domain and light source locations, (a) Illus- 
trates a layout of the inverse problem setting in 2-D. The circular disk Q corresponds to a horizontal 
cross-section of the hemisphere part of phantom (the supposed "mouse head" in animal experiments), 
The computational domain Oi is a rectangle containing f2 inside, 6 light sources are located outside 
of the computational domain. Because of limitations of our device, we use only three locations of 
the light source along one line (number 1,2,3) to model the source xq running along the straight line 
L. Light sources numbered 1,4, 5, and 6 are used to construct an approximation for the tail function, 
(b) Depicts the computational domain Oi = {(a:,2:),|x| < 5.83, |2:| < 5.83} (unit:mm) and its rect- 
angular meshes for tail functions and inverse calculations for the numerical method of this paper. 
Actual mesh is much dense than these displayed. The diagram is not scaled to actual sizes. 



latter was sufficient for our goal of imaging of abnormalities. Sources 1,4,5,6 (Fig. 
1(a) ) were used to approximate the tail function as described below. We construct 
that approximation in the domain Vti depicted on Fig. 1(b) This domain is (units 
in millimeters) 

ill := {x = (x, z) : xi = —5.83 < x < X2 = 5.83, zi = —5.83 < z < Z2 = 5.83} . 

(43) 

5.2.1. The first stage of approximating the tail 

Since this stage essentially relies on positions of sources 1,4,5,6, it makes sense to 
describe this stage in detail here, although it was also described in IT]]. The light 
sources are placed as far as possible from the domain r^i, so that the asymptotic ap- 
proximations ()45p . (146 p of solutions are applicable. Let s*^-^^ be the two-dimensional 
vector characterizing the position of the light source number j. First, consider the 
light source number 1. Denote 



^(1) :=5W(x,z,s(i)) =|(x,z)-s« . 

Using (j5]) , Theorem 12.11 and Lemma 12.3^ we obtain for the function u; = In ti 

1 



(44) 



w \ x,z,s 



(1) 



-kS^^^ - In i 2^/2^") - ^ In S^^^ +Poo (x, z) + 0. 



oo. 



(45) 



Since by (03]) w [x, zi, s*^^)) = 99 (x, s^-*^)) ,x G dQ.r\{z = zi\ , we use (jl5]) to ap- 
proximate the unknown function poo {x, zi) as 



Poo {x,zi) = w (^x,zi,s^^^^ + kS [x,zi,s^^^^ + ]^\D.i^ 



IT 



2S (x,zi,s(i)) / ■ 



(46) 



Formula (|i6]) gives the value of poo {x,zi) only at z = zi. Since Qi is a square. 
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we set the first guess for the tail as the one which is obtained from (I46p by simply 
extending the values at z = zi to the entire domain of Oi, 



w (x,z, s^^'^ 



-kS^^^ (x,2,s(i)) - In (2V2^^ -^In^W (x,z,s^^'^^ +Poo {x,zi), 

where the function S [x, z,s^^^^ is given in (j44p . Next, we compute the function 
u (x,z,s(^)) = exp [w [x,z,s^^'))) and get a^^^ {x,z) , {x,z) G Q via (HH), (|i2]) . 

For light sources 4-6, we repeat the above procedure to get a^^^ (x, z) , a^^^ (x, z) 
and a^^^ (x, z) respectively. Then we consider the average coefficient a(x, z) and set 
(see ^) 

a (x, z) := max(a(x, z), /c^). (47) 

Next, we solve the forward problem ([T|), ([2]) for the light s*^^) with this coefficient 
a (x, z) again to get u (x, z, s) , s = |s('^)| . The final approximate tail function ob- 
tained on the first stage is 

, ^ ^'^'uix.Z.'s) , 

Ti{x,z) = ^f^. 48 

5.2.2. The second stage for the tail 

The second stage involves an iterative process that enhances the first approxima- 
tion for the tail (j48p . We describe this stage only briefly here referring to [17] for 
details. In this case we use only one source number 3, s^^^. Recall that js^'^^j = s. 
Let the function a{x,z) be the one calculated in ()17|) . Denote oi (x) := a(x,2;) . 
Next, we solve the following boundary value problem 

Aui — ai (x) ui = 0, x G ill, 

""ilef^i = V5(x,s) ,x G dVti. 

Then, we iteratively solve the following boundary value problems 

/\Wm - am (x) Wm = [am (x) - Om-l (x)] Um-1, m > 2, 
WmldQi = 0. 

We set Um '■= Um-i + Wm- Next, using (|^T]) and (|12|) with u := Um, we find the 
function Om+i (x) . We have computationally observed that this iterative process 
provides a convergent sequence {om (x)} in L2 (ili) . We stop iterations at m := mi, 
where mi is defined via 

l|a„.,-a^,_i||^^(^^) (49) 



ll«mi-l|lL2(ni) 

where e > is a small number of our choice, and the norm in L2 {^1) is understood 
in the discrete sense. Next, assuming that Umi > 0, we set for the tail 

r(x) = ^Inii^, (x). (50) 

Then, using ([50|) . we proceed with calculating of functions g„ as described above. 
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6. Real Data 



We now describe our experimental setup for collecting the optical tomography data 
from an optical phantom. This phantom is a man-made subject that has the same 
optical property as small animals. Such a phantom is a well-accepted standard to 
test reconstruction methods for real applications before animal experiments. 




Figure 2. A photograph of the experimental setup. An optical phanto m is connected with 4 laser fibers, 
the one on the right hand side is movable to locations 1,2, and 3 (Fig. [Tfajp . The gelatin made phantom 
has the shape of the rectangular box with a hemisphere on its top surface. The hemisphere mimics the 
head of a mouse with a mask exposing the crown part of the mouse. A CCD camera mounted above the 
phantom (not shown) provides light intensity measurements of the top surface of the phantom. 

Fig. [2] is a photograph of our measurement setup. The center of the picture is 
the phantom (rectangular box with a hemisphere on top surface) which, in partic- 
ular, contains a hidden inclusion inside (not visible in this photo). The hemisphere 
mimics the mouse head in animal experiments of stroke studies, and hidden inclu- 
sions of ink-mix mimic blood clots. The four needles are laser fibers that provide 
light sources in our experiments. The fiber on the right hand side of Fig. [2] can be 



moved to three other positions which are source positions 1,2,3 on Fig. 1(a) Three 



other fibers represent positions 4,5,6 of the source on Fig. 1(a) A CCD camera is 
mounted directly above the setup (not shown in photo), and the camera focus is 
on the top surface of the phantom. CCD stands for Charge Coupled Device. CCD 
camera with its sensitivity and response range is commonly used for near infrared 
laser imaging of animals. 

The purpose of this experimental setup is to study the feasibility of using our 
numerical method in an animal stroke model. In the case of a real animal, the 
top hemisphere (meshed shape in Fig. 6.2a) is to be replaced by a mouse head 
and the rectangular block of phantom to be replaced by an optical mask filled 
with a matching" fluid (a tissue-like solution with the optical properties similar 
to the animal skin/skull). The image reconstruction should provide the spatial 
distribution of the optical coefficient a(x) (directly related to the blood content) 
in a 2-D cross-section of the animal brain. 

The goal to is to accomplish a noninvasive imaging means to monitor and in- 
vestigate hemodynamic dysfunction during ischemic stroke in animal models. The 
current experimental setup works only for small animals, such as rats, but not for 
human brains due to the limited penetration depth of NIR light through a larger 
size and thicker bone structure of the skull. The methodology developed here, 
however, may be applicable to a neonate head, while further studies are needed to 
confirm such an expectation. 

Our inverse reconstruction is performed in the 2-D plane depicted in Fig. |3(b)| 
with the optical parameter distribution of the medium inside the circle (at the cen- 
ter ofFig. |3(b)D as an unknown coefficient in the photon diffusion model. Should we 
need a diagnosis of a different cross-section of mouse brain in animal experiments. 
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(a) (b) 

Figure 3. The schematic diagrams of our data acquisition process, (a) schematically depicts a 3-D 
phantom and its hidden inclusion (shown in meshed surface) 5mm below the phantom. The light 
source is placed at various locations on the top surface of the rectangular block (see Fig. \ l(a)^ and 
light intensity measurements are taken on the top surface of phantom, (b ) The measurement surface 
by a CCD camera. The data are collected from the surface of the hemisphere as well as from the 
un-shaded area of top rectangle, both area are lif ted fr om the phantom in this drawing for a better 
illustrative purpose. The rectangular figure in Fig . \3 (b'^ also illustrates the middle a 2-D cross-section 
(meshed circle) of the presumed "animal head", at the boundary of which light intensity data are 
collected for the reconstruction. The light intensity at the 2-D cross-section (except for the boundary) 
is obstructed by the top surface of the hemisphere. Light sources are also located in the same plane 
as this cross-section area. The 2-D inverse problem is solved in this cross-section by ignoring the 
dependence on the orthogonal coordinate. 



we need a different optical mask filled with matching fluid, and repeat both the 
experiment and the reconstruction. 




Figure 4. Dimensions of the ph antom , shown in a vertical cross-section at the center of th e pha ntom 
from side view. The circle in Fig. \l(a)\ corresponds to the boundary of the meshed circle in Fig. \3(b )\ This 
also corresponds to the boundary of the 13 mm diameter circle at a top view above the rectangular box (but 
not related to circle or semicircle in this graph). The 5 mm diameter hollow sphere is for the placement of 
inclusions, filled with liquids of different compositions to emulate strokes. The hidden inclusion is 5 mm 
below the top surface. 



The geometry of the phantom (shown in a vertical central cross-section) is de- 
picted in Fig. [4] with dimensions specified. The phantom is shaped by a hemisphere 
(diameter 13 mm) on the top of a cube of 30mm x 30mm x 30mm. A spherical 
hollow of 5 mm diameter is located inside the phantom, with its top 5mm below 
the top surface (Fig. [4]). The location of this hollow is symmetric here but can be 
in any place when needed. One and two spherical hollows of 3mm diameter inside 
the phantom are also used in experiments (see Section 8 for detail). We fill the hoi- 
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low with liquids made of different kinds of ink/intralipid mix to model strokes by 
blood clots. Intralipid is a product for fat emulsion, which mimics the response of 
human or animal tissue to light at wavelengths in the red and infrared ranges. The 
phantom is made of gelatin mixed with the intralipid. The percentage of intralipid 
content is adjustable. So that the phantom has the same optical parameters as the 
background medium of the target animal model. At the location of the inclusion 
(the hollow), we inject ink/intralipid mixed fluids whose optical absorption coef- 
ficients are 2 times, 3 times and 4 times higher than in the background. Also, we 
use the pure black ink to test our reconstruction method for the case of the infi- 
nite absorption. Different levels of the inclusion/background absorption ratio are 
used to validate our method for its ability for different blood clots. Light intensity 
measurements are taken directly above the semicircle in Fig. HI We note that Fig. 



1(a) shows a top view of the layout but Fig. S] shows a cross-section from side view. 
The disk region in Fig. 1(a) corresponds to the meshed area in Fig. |3(b)| as weh 



as a horizontal cross-section of the 13 mm hemisphere. The reason we use such a 
geometry is that we need to test our capability to reconstruct inclusions at several 
different depth (moving up and down vertically as in Fig. U]) as well as at different 
locations (moving horizontally in Fig. |3(b)] ). 



7. Processing Real Data 

As it is typical when working with imaging from real data, we need to make several 
steps of data pre-processing before applying our inverse algorithm. 



7.1. Computational domains 

We need to use several computational domains. The computational domains 



and meshes used in our numerical calculation involve four domains. Fig. 5(a) 



Fig. 5(c) and Fig. 1(b) are the finite element meshes on each of these domains 
fi, r^O) ^o\^> ^1 respectively. In actual computations we use more refined meshes 
for each domain than those illustrated. Therefore, figures are not scaled to actual 
sizes. 




(a) C (b) (c) Oo\n 

Figure 5. (a) The domain Q = {{x , z) , \/x^~+~z^ < 4.63} (units: mm) with a triangular mesh. The 
domain represents a cross-section of our phantom, and the inclusion inside fl is to be reconstructed, 
(b) The domain Hq = {(1,2), |x| < 23.32, |2| < 23.32} (units: mm). The domain Oq co ntains 
n and light sources. Meshed small square is the domain Hi which is the same as on Fig. \l(b)\ 
Computations of the inverse problem are performed in (c) Shows S7o\^ ihe domain used for 
data processing. To smooth out the measurement noise, equation {Ip with the Dirichlet boundary 
conditions u \gci= ip{x,xo), u \gng= is solved in Qo\Q. Then the smoothed data along dQi is 
used for the inverse problem. 



Here the disk il. is the domain of interest that contains the inclusion, corre- 
sponding to the meshed area in Fig. 3(b) It reflects a cross-section of the hemi- 
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sphere of phantom in Fig. |3(b)[ The hght intensity data used in our computations 
are originated by the measurements at d^l by taking pictures with CCD camera 
on the top of phantom as discussed, shown in Fig. |3(b)[ The CCD camera col- 
lected data is acquired from the surface of the hemisphere and the un-shaded area 
of the top rectangle, and we extract only the needed data for i.e., the cir- 
cle in Fig. |3(b)| by getting the data at these locations as boundary values. Here 

= z) \ \^x'^ + z^ < 4.63mm|. 

On Fig. [5] ilo is a large domain, which can be interpreted as a truncated plane 
[R^ Indeed, one cannot practically solve the forward problem ([l]), ([2]) in the infinite 
plane R^. Light sources are located in f^oX^^- The background simulation and 
calibration of background parameters are performed in f^o- To smooth out the 
noise in the data, equation ([T]) is solved in Qo\^ for each source position, which 
is similar with [24]. Because of (l3|), we use a (x) := k'^ for x G i7o\^- In this 
procedure the Dirichlet boundary condition at dil. is taken from the real data, and 
we use the zero Dirichlet boundary condition at dflQ. As a result, we obtain the 
smoothed Dirichlet boundary condition at dQi for each source location. We solve 
the inverse problem in J^i, and C 



7.2. Data pre-processing and approximations of boundary conditions 

As it was stated in subsection 7.1, to smooth out our measurement data, we solve 
the Dirichlet boundary problem for equation ([T]) with a (x) = /c^ for x G JIqX^ 



(Fig. 5(c) ) for each source location. Namely, 



Au - k'^u = -5 (^x - s« j , X G no\n, (51) 



Here k"^ = 'ifJ-'sfJ-a is the background value and if (x,s*^*^) is the experimentally 
measured data for the light source s*^*) . Let ^ (x, s^*-*) be the trace of the solution of 
the Dirichlet boundary value problem (|5ip at x G dQi. We solve the inverse problem 
in the square ili with the smoothed data u \dni= ^(x, s*^*)) . We have shown 



m 



23] that solving inverse problems in the original domain 17 and extended 



domain are equivalent mathematically. But numerically the noisy component 
in ^(x, s^*^) is much smaller than in if (x, s*^*)) . This smoothing effect takes place 
because the inverse problem is solved in ili. 



7.3. The forward problem and calibration 

We now address the question on which value of k'^ one should use when working with 
the real data. The optical properties of the background of the phantom (without 
the hidden inclusion) are known theoretically from the concentration of intralipid in 
the mix. However, there is a discrepancy between the theoretical value and actual 
measurements. Before we solve the inverse problem to image hidden inclusions, 
we calibrate our model by adjusting the background value of k'^ to the real data 
measured for the reference medium, which is the phantom without inclusion, i.e. 
the hollow is filled with the same intralipid solution as that in the phantom itself. 
First, we numerically solve the forward problem with the source position s^^^ in 
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the domain f^o without any inclusion, 

An- = -y4(5 (^x- s(^)) ,x G J^o, (52) 

u IdUo = 0, 

where k"^ = ?>^'g^a is the background value. Then we calibrate the parameter 
(fixing ^g) as well as the amplitude ^ > of light source in our model (I52p to 
match the measured light intensity for s*^^^ for the uniform background. We do not 
know the number A. Thus, we choose A in such a way that Ucomp (^max?^ ) ~ 
Umeas (xmaxjS^^)) . Here x^ax is the brightest point, i.e. the far right point on dO.^ 
being closest to the light source. Also, Ucomp (xmax, s'^^)) and Umeas (xmaxiS^^^) are 
computed and measured light intensities respectively. Next, we should approximate 
the constant k^. To do this, we take another sampling point Xmin with the minimal 
light intensity, which is the farthest left point on d^l. Next we consider ratios 

Rcomp ) 5 ^measj where 



73 2^ ^comp (Xmax)S^ ^) „ Umeas (^max) 

^cornp j — 7 TTYV' ^' 



' ) "'vaeas y-'^iaim^ ) 

These ratios are independent on the number A in (j52p . We choose k"^ such that 
Rcomp {k'^) ~ Rmeas- As a result, the calibrated value of k'^ was k"^ = 2.403. This 
computed value matches quite well the theoretical value of 2.4 of the intralipid 
solution we have used. 



8. 



Results Of the Reconstruction 



Let 

o-ind — ^ (x) be the value of a (x) inside the inclusion and a}j — k"^ — 2.403 
be the value of the coefficient a (x) in the background, which was computed in 
subsection 7.3. Our ratios aind/o-b where 



O'incl 

at 



2,3,4,00. 



(53) 



The value aind/cib = oo means that the inclusion was filled with a black absorber, 
i.e. black ink. Table [1] gives a detailed information of our real data with different 
positions, numbers, diameters and contrast of inclusions. "Y" means that the ex- 
periment was done for this setting, and "-" means that the experiment was not 
done. Note that in Group 3 we have imaged two different inclusions simultaneously. 



Inclusion Groups 


Positions 


Diameters 


Ratio2 


Ratio3 


Ratio4 


Ratiooo 


Group 1 


@ 


5mm 


Y 


Y 


Y 


Y 


Group 2 




3mm 


Y 


Y 


Y 




Group 3 


O 


3mm 


Y 


Y 


Y 





Table 1. Real data 
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As described in sub-subsection 5.2.1, we construct the "asymptotic tail" on the 
first stage using Ught sources 1,4,5, and 6 as an initial approximation, to be followed 
by other subroutines for further refinements. The image of the function a{x, z) 
in ()47p is depicted on Fig. [U for a phantom where the theoretical value of the 
inclusion/background contrast (f53l) was ainci/at = 3. Images from the asymptotic 
tails for other phantoms listed in Table 8.1 were similar. 




Figure 6. The function a (x, z) in (^7| ) is depicted. This function was obtained on the first stage procedure 
for the tail (sub-subsection 5.2.1) as an initial approximation. The initial reconstruction is obtained from 
a phantom where the theoretical value of the inclusion/background contrast (5.3V was ai„cl/o.b = 3. 

Figures El- Fig. [9] depict the reconstructed images from our real data. Fig. [7| 
depicts the 3D plots of reconstructed images of group 1 for contrast values of: 2:1, 
3:1, 4:1 and oo : 1, from the left to the right respectively. Fig. [8]- Fig. [9] show 2:1, 
3:1, 4:1, from left to right respectively for groups 2,3. In all our examples e = 
in (j49p . The finally reconstructed results for contrasts a^^maxa(x) are listed in 
Table [2j Note that our above reconstruction algorithm does not use any knowledge 
of neither the location of the inclusion, nor the contrast value aind/ab- 




Figure 7. Reconstructed inclusion contrast of data group 1, actual contrasts are 2:1, 3:1, 4:1 and oo : 1, 
from left to right respectively. The transparent frames show the theoretical values of inclusion /background 
contrasts of actual inclusions, which are made of different ink-intralipid mix. 




Figure 8. Reconstructed inclusion contrast of data group 2, actual contrasts are 2:1, 3:1, 4:1, from left to 
right respectively. The transparent frames show the theoretical values of inclusion/background contrasts 
of actual inclusions, which are made of different ink-intralipid mix. 
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Figure 9. Reconstructed inclusion contrast of data group 3, actual contrasts are 2:1, 3:1, 4:1, from left to 
right respectively. The transparent frames show the theoretical values of inclusion/background contrasts 
of actual inclusions, which are made of different ink-intralipid mix. 





The true contrast aind/db 


^ max a (x) 


Relative Error 


Group 1 


2 


2.11 


0.056 


3 


2.9 


0.032 


4 


4.22 


0.057 


oo 


6.69 


unknown 


Group 2 


2 


2.33 


0.167 


3 


3.29 


0.0986 


4 


4.57 


0.143 


Group 3 


2 


2.29 


0.148 


3 


3.49 


0.164 


4 


4.58 


0.147 



Tabic 2. Reconstructed values of the contrast ^ max a (x) within imaged inclusions and they relative errors. 
Of, — 2.403, compare with (8.1). 



9. Summary 

We have worked with a set of real data, using the approximately globally convergent 
numerical method of [17] for a Coefficient Inverse Problem for an elliptic equation. 
These data mimic imaging of clots in the heads of a mouse. We have introduced a 
new concept of the approximate global convergence property and have established 
this property for the discrete case of a finite number of finite elements, unlike the 
continuous case of [ITj]. Figures d- Fig. [9] as well as Table [2] demonstrate that our 
are quite accurate, including both inclusion/background contrasts and locations of 
inclusions. Note that accurate values of contrasts are usually hard to reconstruct 
via locally convergent algorithms. On the other hand, these values are especially 
important for our target application, since they might be used for monitoring stroke 
treatments. Therefore, following the concept of Steps 1-6 of subsection 4.1, we 
conclude that our approximate mathematical model of subsection 4.3 is a valid 
one. 
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